SEIQRHF Network Model
Luis Chaves
## Warning: package 'tibble' was built under R version 3.6.2## Warning: package 'purrr' was built under R version 3.6.2## Warning: package 'lubridate' was built under R version 3.6.2## Warning: package 'broom' was built under R version 3.6.2## Warning: package 'gt' was built under R version 3.6.2## Warning: package 'devtools' was built under R version 3.6.2Sourcing custom modules
source("../Scripts/SEIQHRFNetModules.R")Setting very basic network
n = 200
nw = network.initialize(n = n, directed = FALSE)Getting nodal attributes from data
ageData = read.csv("../Data/age_and_sex.csv") %>% rename(Age = V1, Gender = V2, ID = X)
# I'm just taking Camp 1 data as it seems more complete
campParams = read.csv("../Data/camp_params.csv") %>% rename(Pop_Proportion = Value) %>%
filter(Camp == "Camp_1", Variable == "Population_structure")
campParams$Age = gdata::drop.levels(campParams$Age)
ageGroups = campParams %>%
select(Age) %>% as.matrix()
# I'm assuming age groups to be left inclusive
# and right exclusive but probably does not matter tooo much
ageData$ageGroup = cut(ageData$Age, breaks = c(0,10,20,30,40,50,60,70,80, Inf))
plot(campParams$Age, campParams$Hosp_given_symptomatic)
plot(campParams$Age, campParams$Critical_given_hospitalised)
plot(density(ageData$Age))
paramsFromData = list()
paramsFromData$age.dist = ageData$ageGroup
# the two below are in same order as age groups
paramsFromData$rates.byAge = data.frame(AgeGroup = levels(ageData$ageGroup),
hosp.rate = campParams$Hosp_given_symptomatic,
fat.rate = campParams$Critical_given_hospitalised)Setting network structure based on how refugees are allocated to tents
Based on Tucker model (from Manchester U.) Each individual is a member of a household that occupies either an isoboxor a tent. Isoboxes are prefabricated housing units with a mean occupancyof 10 individuals. Tents have a mean occupancy of 4 individuals. A total of 8100 individuals occupy isoboxes, and 10,600 individuals occupy tents. The exact occupancy of each isobox or tent is drawn from a Poisson distribution, and individuals are assigned to isoboxes or tents randomly without regard to sex or age. This is appropriate because many people arrive at Moria travelling alone, and thus isoboxes or tents may not represent family units
prop.isobox = round(8100/(8100+10600),2)
prop.tent = 1 - prop.isobox
stopifnot(round(n*prop.isobox+n*prop.tent)==n)
# residence = c(rep("isobox", n*prop.isobox), rep("tent", n*prop.tent))Could also think of this problem as people in tents with tent IDs
# how many tents are there
tent.capacity = 4
iso.capacity = 10
num_of_tents = round(n*prop.tent/tent.capacity)
num_of_iso = round(n*prop.isobox/iso.capacity)
housing = c(rep(paste0("tent", 1:num_of_tents),tent.capacity),
rep(paste0("isobox", 1:num_of_iso), iso.capacity))
residence = c(rep("tent",num_of_tents*tent.capacity),
rep("isobox", num_of_iso*iso.capacity))
## If NA (unallocated) put in tent)
if (length(housing)<n){
remaining = n-length(housing)
print(paste(remaining, "nodes were left unassigned. Assigning them to new tent."))
housing[length(housing)+(1:remaining)] = paste0("tent", num_of_tents+1)
residence[length(residence)+(1:remaining)] = paste0("tent")
num_of_tents = num_of_tents+1
}
table(housing)## housing
## isobox1 isobox2 isobox3 isobox4 isobox5 isobox6 isobox7 isobox8 isobox9 tent1
## 10 10 10 10 10 10 10 10 10 4
## tent10 tent11 tent12 tent13 tent14 tent15 tent16 tent17 tent18 tent19
## 4 4 4 4 4 4 4 4 4 4
## tent2 tent20 tent21 tent22 tent23 tent24 tent25 tent26 tent27 tent28
## 4 4 4 4 4 4 4 4 4 4
## tent29 tent3 tent4 tent5 tent6 tent7 tent8 tent9
## 4 4 4 4 4 4 4 4table(residence)## residence
## isobox tent
## 90 116Set vertex attribute to housing or residence
nw = set.vertex.attribute(nw, "housing", housing[1:n])
nw = set.vertex.attribute(nw, "residence", residence[1:n])Settting nodal attributes
nw = set.vertex.attribute(nw, "age", sample(as.vector(paramsFromData$age.dist),n))Explanation of the formation terms (documentation can be found by running help(edge.terms) and choosing the ergm option). We’ll explain here what some basic terms are and should be:
edges: This term adds one network statistic equal to the number of edges (i.e. nonzero values) in the network. For undirected networks, edges is equal to kstar(1); for directed networks, edges is equal to both ostar(1) and istar(1).concurrent: This term adds one network statistic to the model, equal to the number of nodes in the network with degree 2 or higher. The optional term attrname is a character string giving the name of an attribute in the network’s vertex attribute list. If this is specified then the count is the number of nodes with ties to at least 2 other nodes with the same value for that attribute as the index node. This term can only be used with undirected networks. *isolates: This term adds one statistic to the model equal to the number of isolates in the network. For an undirected network, an isolate is defined to be any node with degree zero. For a directed network, an isolate is any node with both in-degree and out-degree equal to zero.meandeg— Mean vertex degree: This term adds one network statistic to the model equal to the average degree of the vertices. Note that this term is a constant multiple of both edges and density.degree(d, attrname)— Degree: The d argument is a vector of distinct integers. This term adds one network statistic to the model for each element in d; the ith such statistic equals the number of nodes in the network of degree d[i], i.e. with exactly d[i] edges. The term attrname is a character string giving the name of an attribute in the network’s vertex attribute list. If this is specified then the degree count is the number of nodes with the same value of the attribute as the ego node. This term can only be used with undirected networks.nodemix(attrname, base = NULL)— Nodal Attribute Mixing: The attrname ar- gument is a character string giving the name of a categorical attribute in the network’s vertex attribute list. This term adds one network statistic to the model for each possible pairing of attribute values.The statistic equals the number of edges in the network in which the nodes have that pairing of values. In other words, this term produces one statistic for every entry in the mixing matrix for the attribute. The ordering of the attribute values is alphabetical (for nominal categories) or numerical (for ordered categories). The optional base argument is a vector of integers corresponding to the pairings that should not be included. If base contains only negative integers, then these integers correspond to the only pairings that should be included. By default (i.e., with base = NULL or base = 0), all pairings are included.nodematch(attrname, diff = FALSE, keep = NULL)— Uniform homophily and differential homophily: The attrname argument is a character string giving the name of an attribute in the network’s vertex attribute list. When diff = FALSE, this term adds one network statistic to the model, which counts the number of edges (i, j) for which attrname(i) == attrname(j). When diff = TRUE, p network statistics are added to the model, where p is the number of unique values of the attrname attribute. The kth such statistic counts the number of edges (i, j) for which attrname(i) == attrname(j) == value (k), where value(k) is the kth smallest unique value of the attribute. If set to non-NULL, the optional keep argument should be a vector of integers giving the values of k that should be considered for matches; other values are ignored (this works for both diff = FALSE and diff = TRUE. For instance, to add two statistics, counting the matches for just the 2nd and 4th categories, use nodematch with diff = TRUE and keep = c(2,4).density:This term adds one network statistic equal to the density of the network. For undirected networks, density equals kstar(1) or edges divided by n(n-1)/2; for directed networks, density equals edges or istar(1) or ostar(1) divided by n(n-1).nodefactor(attrname, base = 1)— Main effect of a factor attribute: The attrname argument is a character string giving the name of a categorical attribute in the network’s vertex attribute list. This term adds multiple network statistics to the model, one for each of (a subset of ) the unique values of the attrname attribute. Each of these statistics gives the number of times a vertex with that attribute appears in an edge in the network. In particular, for edges whose endpoints both have the same attribute value, this value is counted twice. To include all attribute values is usually not a good idea, because the sum of all such statistics equals twice the number of edges and hence a linear dependency would arise in any model also including edges. Thus, the base argument tells which value(s) (numbered in order according to the sort func- tion) should be omitted. The default value, one, means that the smallest (i.e., first in sorted order) attribute value is omitted. For example, if the “fruit” factor has levels “orange”, “apple”, “banana”, and “pear”, then to add just two terms, one for “apple” and one for “pear”, set “banana” and “orange” to the base (remember to sort the values first) by using nodefactor(“fruit”, base = 2:3). For an analogous term for quantitative vertex attributes, see nodecov.nodecov(attrname)— Main effect of a covariate: The attrname argument is a character string giving the name of a quantitative (not categorical) attribute in the network’s vertex attribute list. This term adds a single network statistic to the model equaling the sum of attrname(i) and attrname(j) for all edges (i, j) in the network. For categorical attributes, see nodesociality- Undirected degree: This term adds one net- work statistic for each node equal to the number of ties of that node. The optional attrname argument is a character string giving the name of an attribute in the net- work’s vertex attribute list that takes categorical values. If provided, this term only counts ties between nodes with the same value of the attribute (an actor-specific ver- sion of the nodematch term). This term can only be used with undirected networks. For directed networks, see sender and receiver. By default, base = 1 means that the statistic for the first node will be omitted, but this argument may be changed to control which statistics are included just as for the sender and receiver terms.
Formation
Formation by residence/housing
formationType = "housing"
if (formationType == "residence"){
formation <- ~edges+
# concurrent+
# nodefactor("residence")+ # tent stat
nodematch("residence") # amount of interaction with same class nodes
mean_degree.iso = 10
mean_degree.tent = 4
residence.iso = sum(residence == "isobox")*mean_degree.iso
residence.tent = sum(residence == "tent")*mean_degree.tent
mean_degree = (sum(residence == "isobox")*mean_degree.iso+sum(residence == "tent")*mean_degree.tent)/length(residence)
concurrent_percentage = 0.1 # % of nodes (people) with a degree of 2 or larger
edges = n*mean_degree/2#number of expected edges
# concurrent_nodes = n*concurrent_percentage
residence.match = edges*.8 # 80% of connection people of same class
target.stats = c(edges,
# concurrent_nodes,
# residence.iso,
# residence.tent,
residence.match)
} else {
formation <- ~edges+
# nodefactor("residence")+ # iso stat
nodematch("housing", diff = FALSE) # amount of interaction with same class nodes
external.friends = 0
#' In the eq. below the first RH term is the mean degree per isobox/tent
mean_degree.iso = mean(table(housing)[grepl("isobox", names(table(housing)))]) + external.friends # 10
mean_degree.tent = mean(table(housing)[grepl("tent", names(table(housing)))]) + external.friends #4
residence.iso = sum(residence == "isobox")*mean_degree.iso
residence.tent = sum(residence == "tent")*mean_degree.tent
mean_degree = (residence.iso+residence.tent)/length(residence)
edges = n*mean_degree/2#number of expected edge
residence.match = edges*0.9 # % of connection people of same class, this value does not do good when btwn 0.8 and 1
target.stats = c(edges,
# residence.iso,
# residence.tent,
residence.match
)
}d.rate = 0.0001
# coef.diss = dissolution_coefs(dissolution = ~offset(edges)+
# offset(nodematch("housing", diff = FALSE)),
# duration = c(2, #duration of average edge
# 30), #duration of edges in your housing(same as assuming
# # you stay with your people for 6 months
# d.rate = d.rate) # this correspond to external deaths
coef.diss = dissolution_coefs(dissolution = ~offset(edges),
duration = 180, #duration of edges in your housing(same as assuming
# you stay with your people for 6 months
d.rate = d.rate) # this correspond to external deaths
coef.diss## Dissolution Coefficients
## =======================
## Dissolution Model: ~offset(edges)
## Target Statistics: 180
## Crude Coefficient: 5.187386
## Mortality/Exit Rate: 1e-04
## Adjusted Coefficient: 5.224048Building network and properly fitting network to stats
From netest documentation (help(netest)) The edges dissolution approximation method is described in Carnegie et al. This approximation requires that the dissolution coefficients are known, that the formation model is being fit to cross-sectional data conditional on those dissolution coefficients, and that the terms in the dissolution model are a subset of those in the formation model. Under certain additional conditions, the formation coefficients of a STERGM model are approximately equal to the coefficients of that same model fit to the observed cross-sectional data as an ERGM, minus the corresponding coefficients in the dissolution model. The approximation thus estimates this ERGM (which is typically much faster than estimating a STERGM) and subtracts the dissolution coefficients.
The conditions under which this approximation best hold are when there are few relational changes from one time step to another; i.e. when either average relational durations are long, or density is low, or both. Conveniently, these are the same conditions under which STERGM estimation is slowest. Note that the same approximation is also used to obtain starting values for the STERGM estimate when the latter is being conducted. The estimation does not allow for calculation of standard errors, p-values, or likelihood for the formation model; thus, this approach is of most use when the main goal of estimation is to drive dynamic network simulations rather than to conduct inference on the formation model. The user is strongly encouraged to examine the behavior of the resulting simulations to confirm that the approximation is adequate for their purposes. For an example, see the vignette for the package tergm.
## Warning: `set_attrs()` is deprecated as of rlang 0.3.0
## This warning is displayed once per session.##
## ==========================
## Summary of model fit
## ==========================
##
## Formula: TARGET_STATS ~ edges + nodematch("housing", diff = FALSE)
## <environment: 0x7fc5572d3e00>
##
## Iterations: 16 out of 20
##
## Monte Carlo MLE Results:
## Estimate Std. Error MCMC % z value Pr(>|z|)
## edges -4.94254 0.08485 0 -58.25 <1e-04 ***
## nodematch.housing 17.50304 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log-likelihood was not estimated for this fit.
## To get deviances, AIC, and/or BIC from fit `object$fit` run
## > object$fit<-logLik(object$fit, add=TRUE)
## to add it to the object or rerun this function with eval.loglik=TRUE.
##
## Dissolution Coefficients
## =======================
## Dissolution Model: ~offset(edges)
## Target Statistics: 180
## Crude Coefficient: 5.187386
## Mortality/Exit Rate: 1e-04
## Adjusted Coefficient: 5.224048Diagnostics
cores = parallel::detectCores()-1
if (formationType == "residence"){
dx = netdx(est1,
nsims = 3,
nsteps = 180, # simulating 6 months
ncores = cores,
nwstats.formula = ~edges+concurrent+nodefactor("residence", levels = T)+
nodematch("residence", diff = T)+nodemix("residence"))
} else {
dx = netdx(est1,
nsims = 10,
nsteps = 180, # simulating 6 months
ncores = cores,
nwstats.formula = ~edges+nodefactor("residence", levels = T)+
nodematch("housing"))
}##
## Network Diagnostics
## -----------------------
## - Simulating 10 networks
## - Calculating formation statistics
## - Calculating duration statistics
## - Calculating dissolution statistics
## if (length(dx$stats[[1]])<10){
plot(dx)
par(mfrow = c(1,2))
plot(dx, "duration")
plot(dx, "dissolution")
}
dx## EpiModel Network Diagnostics
## =======================
## Diagnostic Method: Dynamic
## Simulations: 10
## Time Steps per Sim: 180
##
## Formation Diagnostics
## -----------------------
## Target Sim Mean Pct Diff Sim SD
## edges 662.136 634.311 -4.202 14.595
## nodefactor.residence.isobox NA 777.088 NA 19.219
## nodefactor.residence.tent NA 491.533 NA 17.385
## nodematch.housing 595.922 494.102 -17.086 7.675
##
## Dissolution Diagnostics
## -----------------------
## Target Sim Mean Pct Diff Sim SD
## Edge Duration 180.000 49.337 -72.591 37.590
## Pct Edges Diss 0.006 0.006 -0.130 0.003plot(dx)
Running epidemic
param = param.net(act.rate.se = 10,
inf.prob.se = 0.02,
act.rate.si = 10,
inf.prob.si = 0.05,
act.rate.sq = 2.5,
inf.prob.sq = 0.02,
ei.rate = 1/10,
iq.rate = 1/30, #c(rep(1/30, 60), rep(15/30, 120)), # time varying works
ih.rate = 1/100,
qh.rate = 1/100,
hr.rate = 1/15,
qr.rate = 1/20,
hf.rate = 1/50,
hf.rate.overcap = 1/25,
hosp.cap = 5,
hosp.tcoeff = 0.5,
a.rate = 0,
di.rate = d.rate,
ds.rate = d.rate,
dr.rate = d.rate,
ratesbyAge = paramsFromData$rates.byAge
)
init = init.net(i.num = 3,
r.num = 0,
e.num = 0,
s.num = n - 3,
f.num = 0,
h.num = 0,
q.num = 0
)print(Sys.time()-t0)## Time difference of 54.66581 secsres = as.data.frame(sim1)
# to debug:
# Warning in dat$epi$e.num[at] <- c(0, sum(active == 1 & status == "e")) :
# number of items to replace is not a multiple of replacement length
# The simulation time really goes up with the number of edges## Warning: package 'plotly' was built under R version 3.6.2res %>% select(s.num, e.num, i.num, q.num, h.num, r.num, f.num, num, time) %>%
group_by(time) %>% summarise_all(~mean(.)) %>%
pivot_longer(-time) %>% ggplot(aes(x = time, y = value, color = name))+
geom_line(size = 1)+scale_color_brewer(palette = "Set1")
ggplotly(res %>% select(s.num, e.num, i.num, q.num, h.num, r.num, f.num, num, time) %>%
group_by(time) %>% summarise_all(~mean(.)) %>%
pivot_longer(-time) %>% ggplot(aes(x = time, y = value, color = name))+
geom_line(size = 1)+scale_color_brewer(palette = "Set1"))Plot by age groups
res %>% select(contains("i.num.age"), time) %>% group_by(time) %>% summarise_all(~mean(.)) %>%
pivot_longer(-time) %>% ggplot(aes(x = time, y = value, color = name))+
geom_line(size = 1)+scale_color_brewer(palette = "Set1")
res %>% select(contains("f.num.age"), time) %>% group_by(time) %>% summarise_all(~mean(.)) %>%
pivot_longer(-time) %>% ggplot(aes(x = time, y = value, color = name))+
geom_line(size = 1)+scale_color_brewer(palette = "Set1")
ggplotly(res %>% select(starts_with("num.age"), time) %>% group_by(time) %>% summarise_all(~mean(.)) %>%
pivot_longer(-time) %>% ggplot(aes(x = time, y = value, color = name))+
geom_line(size = 1)+scale_color_brewer(palette = "Set1"))For diagnostics
get_times <- function(simulation.object) {
sim <- simulation.object
for (s in 1:sim$control$nsims) {
if (s == 1) {
times <- sim$times[[paste0("sim", s)]]
times <- times %>% mutate(s = s)
} else {
times <- times %>% bind_rows(sim$times[[paste("sim",
s, sep = "")]] %>% mutate(s = s))
}
}
times <- times %>%
mutate(infTime = ifelse(infTime < 0, -5, infTime),
expTime = ifelse(expTime < 0, -5, expTime)) %>%
mutate(incubation_period = infTime - expTime,
illness_duration = recTime - expTime,
illness_duration_hosp = dischTime - expTime,
hosp_los = dischTime - hospTime,
quarantine_delay = quarTime - infTime,
survival_time = fatTime - infTime) %>%
select(s,
incubation_period,
quarantine_delay,
illness_duration,
illness_duration_hosp,
hosp_los,
survival_time) %>%
pivot_longer(-s, names_to = "period_type", values_to = "duration") %>%
mutate(period_type = factor(period_type,
levels = c("incubation_period",
"quarantine_delay",
"illness_duration",
"illness_duration_hosp",
"hosp_los",
"survival_time"),
labels = c("Incubation period",
"Delay entering isolation",
"Illness duration",
"Illness duration (hosp)",
"Hospital care required duration",
"Survival time of case fatalities"),
ordered = TRUE))
return(times)
}times = get_times(sim1)
times %>% filter(duration <= 30) %>% ggplot(aes(x = duration)) +
geom_density() + facet_wrap(period_type ~ ., scales = "free_y") +
labs(title = "Duration frequency distributions", subtitle = "Baseline simulation")
library(ggnet)
if (n<400){
# some mode options are: circle, kamadakawai, fruchtermanreingold
plot(sim1, "network", at = 15,
vertex.col = "age", legend = T, mode = "fruchtermanreingold",
main = "Age")
if (formationType == "residence") {
plot(sim1, "network", at = 15,
vertex.col = "residence",
legend = T, mode = "kamadakawai", main = "Residence")
} else {
plot(sim1, "network", at = 15,
vertex.col = "housing",
legend = T, mode = "kamadakawai", main = "Housing")
}
}
#try and plot this by housing
nwt = get_network(sim1)
net_at_1 = network.collapse(nwt, at = 1)
library(intergraph)
graph = asIgraph(net_at_1)
if (n<=50){ ## plot of many nice network graphs if the network is small enough
library(igraph)
library(RColorBrewer)
clp = cluster_label_prop(graph)
pal = brewer.pal(length(unique(V(graph)$housing)), "Accent")
# tex2num = 1:length(unique(V(graph)$housing))
layouts = "layout_as_star(), layout_as_tree(), layout_in_circle(), layout_nicely(), layout_on_grid(), layout_randomly(), layout_with_dh(), layout_with_fr(), layout_with_gem(), layout_with_graphopt(), layout_with_kk(), layout_with_lgl(), layout_with_mds()"
layouts = strsplit(layouts, "(),", fixed = T)[[1]]
layouts[length(layouts)] = strsplit(layouts[length(layouts)], "\\()")[[1]]
# V(graph)$color = pal[tex2num]
# I think laayout_with_graphopt is the best
for (lay in layouts){
l = eval(parse(text = paste0(lay, "(graph)")))
print(lay)
plot(
clp,
graph,
vertex.label = NA,
layout = l,
edge.color = "gray50",
legend = T)
}
}Network as a heatmap (with adjacency matrix)
library(igraph)##
## Attaching package: 'igraph'## The following object is masked from 'package:plotly':
##
## groups## The following objects are masked from 'package:network':
##
## %c%, %s%, add.edges, add.vertices, delete.edges, delete.vertices,
## get.edge.attribute, get.edges, get.vertex.attribute, is.bipartite,
## is.directed, list.edge.attributes, list.vertex.attributes,
## set.edge.attribute, set.vertex.attribute## The following objects are masked from 'package:lubridate':
##
## %--%, union## The following objects are masked from 'package:dplyr':
##
## as_data_frame, groups, union## The following objects are masked from 'package:purrr':
##
## compose, simplify## The following object is masked from 'package:tidyr':
##
## crossing## The following object is masked from 'package:tibble':
##
## as_data_frame## The following objects are masked from 'package:stats':
##
## decompose, spectrum## The following object is masked from 'package:base':
##
## unionadj = as_adjacency_matrix(graph, sparse = F)
colnames(adj) = V(graph)$housing
rownames(adj) = V(graph)$housing
library(pheatmap)
adj = adj[order(rownames(adj)), order(colnames(adj))]
pheatmap(adj,
color = c("grey50","black"),
border_color = "white",
angle_col = 45,
angle_row = 45,
fontsize = 6,
legend_breaks = c(0,1),
legend = F,
cluster_rows = F,
cluster_cols = F,
show_rownames = ifelse(n<100, T, F),
show_colnames = ifelse(n<100, T, F))
table(V(graph)$housing)##
## isobox1 isobox2 isobox3 isobox4 isobox5 isobox6 isobox7 isobox8 isobox9 tent1
## 10 10 10 9 9 9 9 9 9 4
## tent10 tent11 tent12 tent13 tent14 tent15 tent16 tent17 tent18 tent19
## 4 4 4 4 4 4 4 4 4 4
## tent2 tent20 tent21 tent22 tent23 tent24 tent25 tent26 tent27 tent28
## 4 4 4 4 4 4 4 4 4 4
## tent29 tent3 tent4 tent5 tent6 tent7 tent8 tent9
## 4 4 4 4 4 4 4 4detach("package:igraph")network GIF
library(intergraph)
library(igraph)##
## Attaching package: 'igraph'## The following object is masked from 'package:plotly':
##
## groups## The following objects are masked from 'package:network':
##
## %c%, %s%, add.edges, add.vertices, delete.edges, delete.vertices,
## get.edge.attribute, get.edges, get.vertex.attribute, is.bipartite,
## is.directed, list.edge.attributes, list.vertex.attributes,
## set.edge.attribute, set.vertex.attribute## The following objects are masked from 'package:lubridate':
##
## %--%, union## The following objects are masked from 'package:dplyr':
##
## as_data_frame, groups, union## The following objects are masked from 'package:purrr':
##
## compose, simplify## The following object is masked from 'package:tidyr':
##
## crossing## The following object is masked from 'package:tibble':
##
## as_data_frame## The following objects are masked from 'package:stats':
##
## decompose, spectrum## The following object is masked from 'package:base':
##
## unionlibrary(animation)
nwt = get_network(sim1)
ani.record(reset = TRUE) # clear history before recording
for (at in c(1:30)){
net_at_1 = network.collapse(nwt, at = at)
graph = asIgraph(net_at_1)
adj = as_adjacency_matrix(graph, sparse = F)
colnames(adj) = V(graph)$housing
rownames(adj) = V(graph)$housing
adj = adj[order(rownames(adj)), order(colnames(adj))]
pheatmap(adj,
color = c("grey50","black"),
border_color = "white",
angle_col = 45,
angle_row = 45,
fontsize = 6,
legend_breaks = c(0,1),
legend = F,
cluster_rows = F,
cluster_cols = F,
show_rownames = ifelse(n<100, T, F),
show_colnames = ifelse(n<100, T, F))
ani.record()
}
oopts = ani.options(interval = 0.5)
ani.replay()
# saveHTML(ani.replay(), img.name = "record_plot")
table(V(graph)$housing)##
## isobox1 isobox2 isobox3 isobox4 isobox5 isobox6 isobox7 isobox8 isobox9 tent1
## 9 9 8 9 9 8 6 5 8 4
## tent10 tent11 tent12 tent13 tent14 tent15 tent16 tent17 tent18 tent19
## 3 4 2 4 3 4 4 3 4 4
## tent2 tent20 tent21 tent22 tent23 tent24 tent25 tent26 tent27 tent28
## 4 3 3 4 4 4 4 4 3 3
## tent29 tent3 tent4 tent5 tent6 tent7 tent8 tent9
## 3 3 4 4 4 4 4 4detach("package:igraph")Network plot by clusters
Ideally would look like 
Get degree distribution
Make animation of network
# dynamic animation of network
if (n <=200 & formationType == "residence"){
library(ndtv)
nw = get_network(sim1)
nw = color_tea(nw, verbose = F)
slice.par <- list(start = 1, end = 25, interval = 1,
aggregate.dur = 1, rule = "any")
render.par <- list(tween.frames = 10, show.time = FALSE)
plot.par <- list(mar = c(0, 0, 0, 0))
compute.animation(nw, slice.par = slice.par, verbose = TRUE)
residence <- get.vertex.attribute(nw, "residence")
residence.shape <- ifelse(residence == "isobox", 4, 50)
residence.color = ifelse(residence == "isobox", "red", "green")
age <- get.vertex.attribute(nw, "age")
render.d3movie(
nw,
render.par = render.par,
plot.par = plot.par,
vertex.sides = residence.shape,
vertex.col = residence.color,
edge.col = "darkgrey",
vertex.border = "lightgrey",
displaylabels = FALSE,
filename = paste0(getwd(), "/movie.html"))
}